from math import *
from numpy import *
import matplotlib.pyplot as plt

e = 1.6021*10**-19
m = 9.1091*10**-31
h = 6.6256*10**-34

d1= 2.13*10**-10
d2=1.23*10**-10 #m

L = 0.135    #m
s_L = 0.001

V = [3000, 3500, 4000, 4500, 5000]
s_V = 100
s_D = 0.5/1000  # 0.5 mm de incertidumbre

medida_3kV = ([28.77, 29.48, 28.03, 27.43], [51.27, 51.70, 51.70, 51.70]) #(D1, D2) en mm
medida_35kV = ([27.87, 28.35, 25.90, 27.43], [47.08, 46.63, 47.32, 47.30])
medida_4kV = ([24.65, 25.61, 24.03, 24.53], [43.20, 44.35, 42.80, 42.47])
medida_45kV = ([24.50, 23.57, 23.67, 23.55], [40.78, 40.55, 41.33, 40.98])
medida_5kV = ([22.45, 23.39, 21.91, 22.22], [39.23, 39.78, 39.33, 38.85])

medida_3kV = (float(mean(medida_3kV[0])/1000), float(mean(medida_3kV[1]))/1000)
medida_35kV = (float(mean(medida_35kV[0])/1000), float(mean(medida_35kV[1]))/1000)
medida_4kV = (float(mean(medida_4kV[0])/1000), float(mean(medida_4kV[1]))/1000)
medida_45kV = (float(mean(medida_45kV[0])/1000), float(mean(medida_45kV[1]))/1000)
medida_5kV = (float(mean(medida_5kV[0])/1000), float(mean(medida_5kV[1]))/1000)

def lambda_asociado(d, D, L):
    return (d*D)/(2*L)

def s_lambda_asociado(d, L, s_D):
    return (d*s_D)/(2*L)

print(f"Lambda_1 (3 kV) = {lambda_asociado(d1, medida_3kV[0], L)} +- {s_lambda_asociado(d1, L, s_D)} m ")
print(f"Lambda_2 (3 kV) = {lambda_asociado(d2, medida_3kV[1], L)} +- {s_lambda_asociado(d2, L, s_D)} m ")
print("------------------------------------------------------------------------")
print(f"Lambda_1 (3.5 kV) = {lambda_asociado(d1, medida_35kV[0], L)} +- {s_lambda_asociado(d1, L, s_D)} m ")
print(f"Lambda_2 (3.5 kV) = {lambda_asociado(d2, medida_35kV[1], L)} +- {s_lambda_asociado(d2, L, s_D)} m ")
print("------------------------------------------------------------------------")
print(f"Lambda_1 (4 kV) = {lambda_asociado(d1, medida_4kV[0], L)} +- {s_lambda_asociado(d1, L, s_D)} m ")
print(f"Lambda_2 (4 kV) = {lambda_asociado(d2, medida_4kV[1], L)} +- {s_lambda_asociado(d2, L, s_D)} m ")
print("------------------------------------------------------------------------")
print(f"Lambda_1 (4.5 kV) = {lambda_asociado(d1, medida_45kV[0], L)} +- {s_lambda_asociado(d1, L, s_D)} m ")
print(f"Lambda_2 (4.5 kV) = {lambda_asociado(d2, medida_45kV[1], L)} +- {s_lambda_asociado(d2, L, s_D)} m ")
print("------------------------------------------------------------------------")
print(f"Lambda_1 (5 kV) = {lambda_asociado(d1, medida_5kV[0], L)} +- {s_lambda_asociado(d1, L, s_D)} m ")
print(f"Lambda_2 (5 kV) = {lambda_asociado(d2, medida_5kV[1], L)} +- {s_lambda_asociado(d2, L, s_D)} m ")
print("========================================================================")

def lambda_teorico(V):
    return h/(float(sqrt(2*m*e*V)))

def s_lambda_teorico(V):
    return (s_V*h)/(2*float(sqrt(2*m*e*V**3)))

print(f"Lambda_teo (3kV) = {lambda_teorico(V[0])} +- {s_lambda_teorico(V[0])}")
print("------------------------------------------------------------------------")
print(f"Lambda_teo (3kV) = {lambda_teorico(V[1])} +- {s_lambda_teorico(V[1])}")
print("------------------------------------------------------------------------")
print(f"Lambda_teo (3kV) = {lambda_teorico(V[2])} +- {s_lambda_teorico(V[2])}")
print("------------------------------------------------------------------------")
print(f"Lambda_teo (3kV) = {lambda_teorico(V[3])} +- {s_lambda_teorico(V[3])}")
print("------------------------------------------------------------------------")
print(f"Lambda_teo (3kV) = {lambda_teorico(V[4])} +- {s_lambda_teorico(V[4])}")
print("========================================================================")

# Hacemos regresion lineal para obtener k con la formula D=k*(1/sqrtV)
D1 = [medida_3kV[0], medida_35kV[0], medida_4kV[0], medida_45kV[0], medida_5kV[0]]
D2 = [medida_3kV[1], medida_35kV[1], medida_4kV[1], medida_45kV[1], medida_5kV[1]]
x = [1/float(sqrt(v)) for v in V]   # x=1/sqrt(V)

# Grafica de D1------------- calculo k1
x = array(x)
D1 = array(D1)

k1 = sum(x * D1) / sum(x**2)
D1_fit = k1 * x
r = corrcoef(x, D1)[0, 1]

plt.scatter(x, D1)
plt.plot(x, D1_fit, color="red")
plt.xlabel("1/sqrt(V)")
plt.ylabel("D1")
plt.title(f"k1 = {k1:.4f}, r = {r:.4f}")
plt.grid()
plt.show()

# Grafica de D2------------- calculo k2
D2 = array(D2)

k2 = sum(x * D2) / sum(x**2)
D2_fit = k2 * x
r = corrcoef(x, D2)[0, 1]

plt.scatter(x, D2)
plt.plot(x, D2_fit, color="red")
plt.xlabel("1/sqrt(V)")
plt.ylabel("D2")
plt.title(f"k2 = {k2:.4f}, r = {r:.4f}")
plt.grid()
plt.show()

# Incertezas de k
s_k1 = 0.0123556
s_k2 = 0.0173088
# Calcular la distancia de red:
def d_red(k):
    return (2*L*h)/(k*float(sqrt(2*m*e)))

def s_d_red(k, s_k):
    return (2*L*h*s_k)/(k**2 *float(sqrt(2*m*e)))

print(f"d1_red = {d_red(k1)} +- {s_d_red(k1, s_k1)} m")
print(f"d2_red = {d_red(k2)} +- {s_d_red(k2, s_k2)} m")
print("========================================================================")

# Calcular constante de planck
def h_exp(d, k):
    return (d*k*float(sqrt(2*m*e)))/(2*L)

def s_h_exp(d, s_k):
    return (d*float(sqrt(2*m*e))*s_k)/(2*L)

print(f"h (d1, k1) = {h_exp(d1, k1)} +- {s_h_exp(d1, s_k1)} J*s")
print(f"h (d2, k2) = {h_exp(d2, k2)} +- {s_h_exp(d2, s_k2)} J*s")